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Abstract. We design consistent discontinuous Galerkin finite element schemes 
for the approximation of the Euler— Korteweg and the Navier-Stokes-Korteweg 
systems. We show that the scheme for the Euler-Korteweg system is energy 
and mass conservative and that the scheme for the Navier-Stokes-Korteweg 
system is mass conservative and monotonically energy dissipative. In this case 
the dissipation is isolated to viscous effects, that is, there is no numerical dis- 
sipation. In this sense the methods is consistent with the energy dissipation 
of the continuous PDE systems. 



1. Introduction 

In this work we propose a new class of finite element methods for the Navier- 
Stokes-Korteweg system which are by design consistent with the energy dissipation 
structure of the problem. The methods are of arbitrary high order of accuracy and 
provide physically relevant approximations free of numerical artifacts. It seems that 
these are the first methods in the literature enjoying these properties. 

Liquid vapour flow occurs in many technical applications and natural phenom- 
ena. A particularly interesting and challenging case is when the fluid undergoes 
phase transition, i.e., there is mass transfer between the phases, which is driven 
by thermodynamics. The applications of these phenomena are extremely varied, 
for example, it is applicable to modelling the fuel injection system in modern car 
engines and also to the study of cloud formation. The modelling of these phenom- 
ena can be traced back to |vdWl iKorOlj . however there remain open questions, for 
example, what is the correct model for the given application at hand. 

The compressible flow of a single substance containing both a liquid and vapour 
phase undergoing a phase transition can be modelled by different techniques. One 
widely used approach for the treatment of these problems, which emerged in the 
last few decades, see |AMW98] and references therein, is the so called diffuse inter- 
face approach. In this philosophy the phases are separated by a (thin) interfacial 
layer across which the fields vary smoothly. The benefit of this approach is that 
there is only one set of PDEs solved on the whole domain whose solution already 
includes the position of the interfacial layer. However, these models must include 
a parameter distinguishing when we are in one phase or another. In most diffuse 
interface models this is a more or less arbitrary indicator function based on the 
mass or volume fraction of one of the constituents. 
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In this contribution we will consider the isothermal Navier-Stokes-Korteweg sys- 
tem which is a diffuse interface model but here the mass density serves as a phase 
indicator, it originates in the work of Korteweg jKorOlj and van der Walls [vdW 
and was derived in modern terminology in DS85, TN92, JLCDOl]. This model 
includes surface tension effects by a third order term in the momentum balance 
which corresponds to a non-local (gradient) contribution in the energy functional. 
Another feature of compressible diffuse interface models is a non-monotone con- 
stitutive relation for the pressure. This corresponds to a non-convex local part of 
the energy, see equation (2.7). As can be seen from Figure [I] the phases of the 




problem (liquid/vapour) are the corresponding regions where the pressure function 
is monotonically increasing. 

The Korteweg type third order term together with the non-monotonicity of the 
pressure function cause several issues in the numerical treatment of this problem. In 
previous numerical studies [JTB02. Die07, BP it has been observed that "classical" 
explicit-in-time finite volume (FV) and discontinuous Galerkin (DG) schemes which 
use standard fluxes used in the computational conservation laws introduce several 
numerical artifacts. 

The first artifact is non-monotonicity of the energy. The Euler-Korteweg model 
is energy conservative over time whereas the Navier-Stokes-Korteweg model is 
monotonically energy dissipative. In the Navier-Stokes-Korteweg model all the 
dissipation is due to viscous effects (see Lemma 2.3). The classical FV and DG 
methods applied to the Navier-Stokes-Korteweg system lead to a non-monotone 
behaviour of the energy. This is mainly due to the fact that these "classical" 
schemes introduce standard diffusion in the mass conservation equation as a stabilis- 
ing mechanism. While for convex energies standard diffusion in fact leads to energy 
dissipation, it may lead to an increase in energy for multiphase flow [Die07( IDGRj . 
In fact standard diffusion is also present in the finite clement method proposed in 
IBP1. 
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The second artifact are so called parasitic currents, i.e., the schemes are not 
well-balanced, as they do not preserve the correct equilibria. Parasitic currents 
occur when equilibrium is approached and the numerical velocity field does not 
vanish uniformly, but in the interfacial layer large velocities whose magnitude is 
dependent on the gridsize and inversely dependent on the width of interfacial layer 
appear |Die07|, §5]. As the interfacial layer is extremely thin this effect cannot be 
neglected in practical computations. 

Both the non-monotone behaviour of the Navier-Stokes-Korteweg energy and 
the parasitic currents are due to numerical regularisation terms which are not 
adapted to the variational structure of the problem, see [PGR] for a study on 
regularisation terms taking into account the underlying variational structure of 
the problem. For previous works on scalar dispersive equations by discontinuous 
Galerkin methods we refer to [CS081 IBTTKlOTl IXSTT] . 

The key target in the work at hand is to consider a high-order DG discretisation 
of the problem which aims at preserving the energy dissipation inequality satisfied 
by the original problem and avoiding the introduction of any artificial diffusion 
terms. By achieving this goal we can treat the case where the system preserves 
the energy exactly. In addition, we can address the case where the system has 
natural dissipation and the energy is diminishing. Our schemes are therefore energy 
consistent in the sense that they are consistent with the energy dissipation structure 
of the Navier-Stokes-Korteweg system. The resulting schemes are free from the 
above mentioned artifacts of other approximating methods in the literature and 
are successful in computing the physically relevant solution. It is to be noted that 
our approach does not hinge on an adaptation of "entropy conservative schemes" 
developed for conservation laws, [Tad03| . The non- monotone pressure function 
makes a direct application of this approach unfeasible in our case. Conservative 
DG schemes for the scalar generalized KdV equation were suggested recently in 
BCKX11 . To achieve our goals we follow a constructive step-by-step approach. 
Motivated by the proof of energy conservation at the continuous level we introduce 
a new mixed formulation for the Navier-Stokes-Korteweg system. This mixed 
formulation will be the basis of our discrete schemes. We first discretise in space by 
employing a DG approach with generic discrete fluxes. Then we specifically identify 
the properties and thus the fluxes which yield energy consistent schemes. Then 
we consider Crank-Nicolson type time discretisation and identify the precise time 
discrete method which is energy consistent. By combining the ideas of space and 
time discretization we obtain the fully discrete schemes with the desired properties. 

The structure of the paper is as follows: In Sj2]we introduce the Navier-Stokes- 
Korteweg model problem as well as some of its conservative properties. We give the 
mixed formulation and necessary notation which will be used throughout the paper. 
In addition we describe the mixed formulation in the broken Sobolev framework 
necessary for the construction of the DG scheme. In ^3] we detail the construction 
of the energy consistent DG scheme initially in the spatially semidiscrete case. We 
then move on to the temporal semidiscrete case in ffl] and combine the results to 
obtain an energy consistent fully discrete scheme in q5j In f|6] we perform various 
numerical experiments to test the convergence, conservativity and computational 
properties of the scheme. 
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2. Model problem, mixed formulation and discretisation 

In this section we formulate the model problem, fix notation and give some basic 
assumptions. Let ft C R d , with d = 1,2, 3 be a bounded domain. We then begin 
by introducing the Sobolev spaces |Cia781 IEva98| 

(2.1) H fe (fi) := {4> € L 2 (0) : D a <f> £ L 2 (0), for |a| < k} , 

which are equipped with norms and semi-norms 

(2-2) IN* := H«llH» ( n) = £ W^Wlm 

\a\<k 

(2-3) and \u\ 2 k := |u£* (n) = ^ IID^H^) 

\ a \=k 

respectively, where a = {ai, ad} is a multi-index, |a| = ^2f—i on and derivatives 

D a are understood in a weak sense. In addition, let 

(2.4) 

Hj := {0 € H x (fi) : 4>\ dn = 0} and H^(fi) := {cf> € [H 1 (fi)] d : (0U) T n = o} 

where n denotes the outward pointing normal to dfl. 

We use the convention that for a multivariate function, u, the quantity Vu is a 
column vector consisting of first order partial derivatives with respect to the spatial 
coordinates. The divergence operator, div , acts on a vector valued multivariate 
function and Am := div (Vu) is the generalised Laplacian operator. We also note 
that when the Laplacian acts on a vector valued multivariate function, it is meant 
componentwise. Moreover, for a vector field v, we denote its Jacobian by Dv. We 
also make use of the following notation for time dependant Sobolev spaces: 



(2.5) L 2 (0,T;H fe (ft)) := Su : [0,T] -> H*(fi) : J ||u(t)||jjdt< 



2.1. Model problem. Consider a fluid in the domain ft with density p and velocity 
v. The Navier-Stokes-Korteweg system is made up of the balances of mass and 
momentum of said fluid, that is, 

d t p + div (pv) = 

(2-6) n / n ,. / \ r—, f \ . ^ A in n x (0,T) 

d t (pv) + div (pv ® v) + Vp{fi) = pAv + 7pVAp 

where p is a non-monotone pressure function, p is a viscosity coefficient and 7 a 
capillarity coefficient. The pressure function p is linked to a double well-potential 
W = W(p) via the relation (Figure llj 

(2.7) P (p) = pW'(p) - W(p). 



Let n be the outward pointing normal to dfl, suppose the system (2.6) is given 
with boundary conditions 

(2.8) v = and (Vp) T n = on dn x (0, T) 
and initial conditions 

(2.9) P(;0)=p°, v(;0)=v°mn 
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for given functions p° G H x (n) and v° e (H 1 ^))"* such that W(p°) € Li(O). The 
system (2.6) conserves mass as well as satisfying a momentum balance together 
with an energy dissipation equality, i.e., 

(2.10) d t (^pda;)=0 

(2.11) d t ( [ pvdx)=p[ (Dv)nds 

(2.12) c 



on 



W{p) 



\Vp\ 2 dxj = -fij \Dv\ 2 dx, 



respectively. The energy dissipation equality is only valid for smooth solutions. In 
case the system permits shocks, they would trigger additional energy dissipation 
and (2.12| would have to be replaced by an inequality. While the first two equalities 
follow by integrating the mass and momentum balance (2.6). The derivation of the 
energy dissipation equality is a little bit more involved. For completeness the result 
is formulated as Lemma 2.3 Moreover, the proof of Lemma 2.3 serves as a guideline 
for the construction of energy consistent discrete schemes. 

2.2. Assumption (finite Helmholtz energy). From hereon in we will assume that 
for a given p we have that W(p) S Li(0, T; Li(Sl)). 

2.3. Lemma. For every smooth solution (p,v) € L 2 (0, T; H 3 (S1)) x L 2 (0, T; H 2 (£l)) d 
of |2li) such that {dtp, d t v) 6 L 2 (0, T; L 2 (fi)) x L 2 (0, T; L 2 (Jl)) d which satisfies the 
boundary conditions (2.8) we have 



(2.13) 



W(p) 



1 I i 2 

nP\ V \ 



l |Vp| 2 d* 



Wvr dx. 



Proof Let us first note that the second equation of (2.6) can be reformulated as 
(2.14) pd t v + div{pv v) - div(pv)v + p\7W'{p) - p,Av - jp\7Ap = 0. 



Multiplying the first equation of (2.6) by W (p) + \ \v\ — 7Ap we see 



(2.15) 



= W'{p)d t p + W'(p) div (pv) + 1 \v\ 2 dtp + i \v\ 2 div (pv) 
— ^Apdtp — 7Apdiv (pv) . 



Then by multiplying (2.14) by v and summing together with (2.15) we obtain 
(2.16) 



= W'(p)d t p + l -d tP \v\ 2 - -/dtpAp + W'ip) div{pv) - l - div{pv) \v\ 2 



7 d\v(pv)Ap + pv J dtV + v 1 div(pv ® v) + pv J WW' (p) 
pv 1 Av — 7/9u T VAp. 



We integrate (2.16) over SI and by Greens formula we have 



(2.17) 



= / W'(p)d tP +\dtp\v\ 2 + pv^dtV + ^d t \V p\ 2 + p\Dv\ 2 dx 



no 



1 



-idtpVp + -pv \v\ — jpvAp + pvW'(p) 



- pv 1 Dvj n ds. 

The boundary integral in (12. 171) vanishes because of the boundary conditions. □ 
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2.4. Remark (stable steady states). The energy dissipation equality gives rise to 



the fact that the (stable) steady states of (2.6) are minimizers of the energy func- 
tional 

(2.18) E\p, v] := J W(p) + \p \v\ 2 + 1 |Vp| 2 dx, 
under the constraint 

(2.19) / pdx = m, 

for some given m > and therefore satisfy the Euler-Lagrange equations 

(2.20) v = 

(2.21) W'(p)-jAp = X, 

where A is the Lagrange multiplier associated with the mass conservation contraint 



(2.191 



Note that (2.21) is equivalent to 

(2.22) = V(W'(p)- 1 Ap) 

= Vp(p) — 7pVAp 

using the relation 

(2.23) Vp(p) = pVW'(p) 
which is readily derived from \2.7\. 



2.5. Classical solvability of the problem. The well-posedness of the Navier- 
Stokes-Korteweg system and similar systems was considered by several authors 
[BGDDJ07) IBDL03) IDD011 IFei02( IHL96) IKot08j . For completeness we will state 
some results. 

2.6. Theorem (existence of a solution to the Euler-Korteweg system [BGDDJ07| ). 
Let s > | + 1 and 

(2.24) H a := H s+1 (R d ) x H s (R d , R d ). 

Suppose the initial data (po,v ) e (p(0),v(0)) + H s where p,v is a special solution 
such that p is bounded away from zero and the Hessian of p, DVp, as we]] as the 
Jacobian of v, Dv are both C([0, T},H s+3 (R d , R dxd )) for some T > 0. Then the 
Euler-Korteweg system admits a unique solution (p,v) <E (p, v) + C 1 ([0, T), H s _ 2 )^ 
C([0,T),H s ) satisfying the initial data (po,v ). 

2.7. Theorem (existence of a solution to the Navier-Stokes-Korteweg system 
[DDOlp . Let B s = i3| 1 (IR d ) denote the homogeneous Besov space. Let p > 
be a reference density such that p'(p) > 0. Suppose also that the initial data po, Vq 
satisfies p - p G B d / 2 , p > c> and v G (B d / 2 ~ 1 ) d . 

Then there exists a T > such that the Navier-Stokes-Korteweg system has a 
unique solution (p,v) with initial data po,v such that p — p G C([0, T), B d / 2 ) n 
L 1 ([0,T),B d /*+ 2 ) and v e C([0,T),(B d / 2 ' 1 )) d D L 1 ([0,T),{B d / 2+1 )) d . 



2.8. Remark. Theorems 12.61 and 12.71 motivate us to construct numerical schemes 
which are adapted to the smooth situation. In particular, enforcing the energy 
dissipation equality proven in Lemma |2.3| 
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2.9. Mixed formulation. To mimic the proof of Lemma [273] at the discrete level 
it will be essential to have at our disposal a numerical formulation in which v 
and t, which depend nonlinearly on the original variables p and pv, are permitted 
as test functions. Indeed, this is our main motivation to reformulate the Navier- 
Stokes-Korteweg system (2.6) as a mixed system of PDEs by the introduction of 
two auxilliary variables, r and q, and using the relation of the pressure function 
and the double well potential (2.23). 

The mixed formulation is then to seek (p,v,r,q) such that 

(2.25) d t p + div (pv) = 

1 2 

(2.26) pd t v + div (pv <g) v) — div (pv) v + pVr — -pV \v\ — pAv = 

(2.27) r- W'(p) + 7 div(q) - - |u>| 2 = 

(2.28) q-Vp = 
which is coupled with the boundary conditions 

(2.29) v = and q J n = 0on ffix (0, T) 
and the initial conditions ( |2.9| . 

2.10. Remark (alternate notation). We note that the second and third term on 
the left hand side of (2.26) can be rewritten as 

(2.30) div (pv ®v) — div (pv) v — p (v 1 V) v, 
which is the standard notation in the incompressible scenario. 

2.11. Discretisation. Let J be a conforming, shape regular triangulation of f2, 
namely, 3? is a finite family of sets such that 

(1) K £ S? implies K is an open simplex (segment for d = 1, triangle for d = 2, 
tetrahedron for d = 3), 

(2) for any K , J € 3? we have that K D J is a full subsimplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) [j Ke ?K = n. 

We use the convention where h : O — > R denotes the meshsize function of 3? , i.e., 

(2.31) /i(cc) := max/ijf, 

where Kk is the diameter of an element K. We let $ be the skeleton (set of common 
interfaces) of the triangulation & and say e € if e is on the interior of Q and 
e g 90 if e lies on the boundary dft. 

2.12. Definition (Broken Sobolev spaces, trace spaces). We introduce the broken 
Sobolev space 

(2.32) R k (3?) := {<t> : <f>\ K E B k (K), for each Kefj, 

similarly for Hj(^) and H^(^). 

We also make use of functions defined in these broken spaces restricted to the 
skeleton of the triagulation. This requires an appropriate trace space 

(2.33) T(£):= J] L 2 (dif)c J] RHdK). 
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Let P P (£T) denote the space of piecewise polynomials of degree p over the trian- 
gulation & we then introduce the finite element spaces 



(2.34) 
(2.35) 
(2.36) 



V := DG(5*,p) = F p {^) 
V := VnHj(^) 
V := V^nHi^) 



to be the usual spaces of (discontinuous) piecewise polynomial functions. For sim- 
plicity we will assume that V is constant in time. 

2.13. Definition (jumps and averages). We may define average and jump operators 
over T (S~) for arbitrary scalar, v G T (<?), and vector valued functions, «£T (<%') d . 

(2.37) " : TV) - \f) 



(2.38) 



(T (<£")) -> {U{S)f 

V !-> |(«kl+«|K 2 ) 



(2.39) 



H : (T(O) 



(2.40) 



: (T (<?))" (L 2 (<?)) 

w ^ (v\K 1 ) 1 n Kl + (v\ K2 ) T n K2 



(2.41) 



(T(<nJ<9ft)) d -> (L 2 (<f)) dxd 



where denotes the outward pointing normal to -fQ. Note that on the boundary 
of the domain dQ the jump and average operators are defined as 



(2.42) 
(2.43) 



vn \\v\\ 



an 



mi 



mi 
v} 



;= v 1 n 



:= v. 



m> 



2.14. Elementwise formulation and discrete fluxes. As a next step towards 
the construction of a numerical scheme we give the elementwise variational for- 
mulation to the problem in mixed form ( 2.25[ )- (2.28). It requires to find (p,v) G 
L 2 (0, T; H 1 ^)) x (L 2 (0, T; Hl(^)) d with (d t p, d t v) G L 2 (0, T; L 2 (^)) x (L 2 (0, T; L 2 (^)) d 
and (r,q) G L 2 (0, T; H x (^)) x L 2 (0, T; H^(^)) such that W'O) G L 2 (0, T; L 2 (^)) 
and 
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(2.44) 



/ (dtp + div (pv)) tp dx + / F 1 (p,v,T,q,ip) ds V^eH 1 ^) 
Jn Jg 

J i^pd t v + div (pv <x> v) — div (pv) v + pVr — ^pV |i>| 2 ^ xdx 

+ f F 2 (p,v,T,q, X )ds + pB(v, X ) V X e(Hj(^)) d 
Jg 

f (r - W'(p) + 7 div (q) - l -\v\ 2 ~) £ d* 

+ / F a (p, 17,^9,0 ds VeeH 1 ^) 
/ (9-Vp) T Cda;+ f F 4 (p,v,T,q,Qds VCeH^) 



where 



Fi,F 3 : H 1 ^) x Hl(^) d x H 1 (^") x H^(^) x H x (^) -> L a (<?) 

(2.45) F 2 : H 1 (^") x H^(^) d x H 1 ^) x H* x Ul(,T) d -> L 2 (#) 

F 4 : H 1 ^) x H*(^) d x H 1 ^) x H^(^) x H^(^) -> L 2 (<?) 

are appropriate choices of elementwise fluxes to be chosen in the sequel to suit 
our purposes; the operators div and V are understood henceforth to be defined 
elementwise and B : (Hj(^)) x (Hj(«^)) -> R is a bilinear form, corresponding to 
a weak formulation of the Laplacian. We will also assume that the fluxes Fl , . . . , F4 
only depend on the traces of their arguments and are linear in the test functions. 
We would like to mention that the spaces for the variational formulation are chosen 



such that all integrals in (2.44) are well-defined. We do not claim that there is 



a wcll-posedness analysis for (2.441 with the given spaces; (2.44 1 serves only as a 



basis to define the spatial discrete DG scheme in the next section. 

3. Development of energy consistent numerical methods - the 

spatially discrete case 

In this section wc will detail the methodology behind the construction of the 
energy consistent finite element scheme. We present our main results which show 



the conditions a generic scheme applied to the variational formulation (2.44) with 
no diffusion (i.e., p, = 0) is mass and energy conservative. If a scheme conserves 
mass for p — this does not change for p 7^ 0, as the mass conservation equation 
is not affected by a reasonable discretization of the viscosity. Moreover if a scheme 
conserves energy for p = 0, for p =/= all energy dissipation is due to viscosity. 
For simplicity we first detail the calculations for the spatially discrete case, then 
construct a temporally discrete scheme. We also give a condition when a scheme 
falling under our framework can also conserve momentum. 

3.1. Spatially discrete scheme. Throughout the calculations in this section we 
will regularly refer to the following proposition. 

3.2. Proposition (elementwise integration). Let 

(3.1) H div (^) := {p G (L 2 (^)) d : divp G L 2 (^)} . 
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Suppose p € H div (^) and (f> £ H 1 ^) then 

(3 2) /_] / div(p)4>dx= 2~] \~ / p J V(f>dx + / ^>p T n/fdi 



A-e,T 



Jn particular we have pGT and <fi £ T {$), and the following identity holds 
(3.3) 

V / WriK ds = [ [p] { ds + / [0] T { p} d,s - / M ds. 

^ JS.R" J^ Jgudn Jguan 



3.3. General numerical scheme. A generic spatially discrete DG formulation 
to the problem in mixed form ( 2.25 )-( 2.28) is to find Phi T h '■ [®,T] — >• V and 
v h : [0, T] -> \ d and q h : [0, T] -> V such that 

(3.4) 

= / (&ph+div(p fc « fc ))*da! + / Fi (p^v^Th^h,^) ds V*eV 

= J {^Phd t Vh + div {p h v h (g) ■u/j) - div (p h v h ) v h + p h \7r h - ^PhV \v h \ 2 ^ Xdx 

+ f F 2 (p h ,v h ,Th,q h ,X) ds + u.B h {v h ,X) VXeV d 
J s 

= J (r h - W'(p h ) + 7 div (g h ) - i |^^| 2 ) Edx 

+ / F 3 (p h ,v h ,T h ,q h ,E) ds VSeV 

0= / (q h ~ Vp h ) J Zdx+ [ F A { Ph ,v h ,r h ,q h ,Z) ds V Z e V, 
Jo J<? 

O , O , 

where : V a x V a -> IR is a discretization of £>. 

3.4. Consistency and conservation. Now we will give abstract properties of 
the fluxes which determine whether the scheme is consistent and conserves mass, 
momentum or energy. 



3.5. Definition (consistency). A generic scheme having the form (3.4) is said to 
be consistent provided 

(3.5) F i (p,v,r,q,-) = 0iori=l,...,A 

for all smooth functions p, t £ L 2 (0, T; H 1 ^)) and v,q £ [L 2 (0, T; H 1 ^))] d . 



3.6. Theorem (conservation). For p — a generic scheme of the form |3.4[ ) con- 
serves: 

(1 ) Mass, that is, 

(3.6) d t 
if and only if 

(3.7) / F 1 {p h ,v h ,r h ,q h ,l)ds = 



p h dx = 



\p h v h \ ds Vp h , T h £ V, v h £ V d , q h £ V, 



where 1 is the the constant element of rl l {3?) which is 1 everywhere. 
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(2) Energy, that is, 
(3.8) d t ( f W( Ph ) + \ Ph \v h \ 2 + I \q h \ 2 dx)=0 



if and only if 

Fi(p h ,v h ,T h ,q h ,T h ) + F 2 (p h ,v h ,T h ,q h ,v h ) + {phT h Vh\ ds = 0, 



s 



(3-9) 

/ F 3 {p hl v h ,T h ,q h ,dtPh)-jDtF 4 (p h ,v h ,T h ,q h ,q h )+~fldtPhqhl ds = 0, 
Jg 

for all p h ,T h : [0,T] -> V,v h : [0,T] -> V d and q h : [0,T] -> V. JVote that 
we use tie notation D t F4(ph,Vh,Th,qh, Z) for the time derivative since 
ph,Vh,Th,qh are time dependent hut Z is independent of time, as V is 
independent of time. 

3.7. Corollary (Energy dissipation). Let be a coercive discretisation of B, then 



for p > a generic scheme of the form {3.4) conserves mass, if and only if (3.7) is 



satished and it satisfies the energy dissipation equality 

(3.10) dt (J Q W(Ph) + \p h \v h \ 2 + \ \q h \ 2 daj) = -pB h (v h ,v h ) < 

if and only if 

/ F 1 (p h ,v h ,T h ,q h ,T h ) + F 2 (p h ,v h ,T h ,q h ,v h ) + lp h T h v h j ds = 0, 

(3.11) 

/ F 3 (p h ,v h ,T h ,q h ,d tPh ) ~ -fD t F 4 (p h ,v h ,T h ,q h ,q h ) +-fld t p h q h j ds = 0, 

holds. 

3.8. Corollary (Consistency, conservation and dissipation). For p = the following 
spatially discrete scheme 

(3.12) 

= / (d tPh + div (p fc » h )) * dx - [ lp h v h j { ds V * e V 

= y {phdtVh + div ® u h ) - div (p h v h )v h + p h \7r h - ^PfeV Iv/J 2 ) Xda; 

- / HfiphX} ds + pB h (v h ,X) VXeV d 
= J (r h - W'{ Ph ) + 7 div(g h ) - i |u h | 2 ) Sdz 

- / )W{S} ds VSeV 

0= / (g /l -V P , l ) T Zd a; + / W T {Z}ds VZeV 

Jf2 J<f 



is consistent and conserves mass (3.6) and energy {3 
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Let Bh(u, w) be the symmetric interior penalty method for the componentwise 
Laplacian given by 
(3.13) 



B h (u,w)= [ Du.Dw dx — f {Dib}:[u] 8 +{ 

Jn Jsuon 



BuMwi 



where : denotes the Frobenius inner product between two d x d matrices, i.e., 
X:Y := trace (X J Y). It is well know for large (enough) a this is a coercive 
discretisation of the componentwise Laplacian. Thus, for p, > 0, tie numerical 



{3.10) 



scheme (3.12) with (3.13) is consistent, conserves mass (3.6) and dissipates energy 



Proof of Theorem |3.6|Lct us first consider the proof of conservation of mass. By 



using \f' = 1 as test function in (3.4)! we want to show 



(3.14) 



= d t ( / p h dx) = / d t phdx = - div(p h v h ) dx - / F 1 (p h .v h ,T h ,q h ,l) ds 

= - {phVh} { 1} dx- Fi(p h ,v hl T h ,q h ,l) ds, 
J g J g 

by Proposition |3.2| with p = phVh and = 1 and noting = on dtt. Hence the 



scheme conserves mass if the condition (3.7 1 is true 



Let us now turn to the conservation of energy. Define 

(3.15) E(p h ,v h ,q h ) := J^W{p h )+ l -p h \v h \ 2 + 1 -\q h \ 2 dx. 
Again we want to show 

(3.16) 0= d t E{p h ,v h ,q h ). 
Explicitly computing the time derivative 

(3.17) 

d t E(p h ,v h ,q h ) = / W'{p h )d t p+ -d t ph \vh\ 2 + Ph(vh) 1 d t v h + j(q h ) J d t q h dx. 



we see 



In view of (3.4)4 and Proposition 3.2 
(3.18) 

d t E(p h ,v h ,q h ) = \ W'(p h )d t p+ -d t ph \vh\ 2 + Ph{vhV d t v h + 7(V<9 t p /l ) T q h dx 



-7 / D t F A (p h ,v h ,T h ,q h ,q h )ds 
Jg 

= I W'(p h )d t p+ T^dtPh \v h \ 2 + Ph(vh) J d t v h - -f(dtPh) 1 divq h dx 

-7 / D t F A (p h ,v h ,T h ,q h ,q h )-l(d t p h )q h j ds, 
Jg 
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as q h J n = on <90. Making use of (3.4)2 an d (3.4) 3 we see 
(3.19) 

d t E(p h , v h , q h ) = I T h d t ph - v h J div(p h v h <8> v h ) + div(p h v h ) \v h \ 2 - puVfJ^Th 
Ju 

+ T^PhVh 1 ^ (\Vh\ 2 ) dx 



Now by (3.4)i and Proposition 3.2 



- / jD t F4,(p h ,v h ,T h ,q h ,q h ) - ^l(d t ph)q h j 
Jg 

- F3(Ph,v h ,T h ,q h ,d t p h ) + F 2 (p h ,v h ,T h ,q h , v h )ds. 
we see, as Vh € \ , 



(3.20) 

d t E(p h ,v h ,q h ) = 



- dxv(p h v h )T h - v h J div(p h v h (g) Vh) + dW(p h v h ) \v h \ 



- PhVhVT h + ^p h v h S7 (\vh\ 



dx 



- / -fD t F 4 (p h ,v h ,T h ,q h ,q h ) - -fl(d t p h )q h } 
Jg 

- F 3 (p h ,v h ,T h ,q h ,d t ph)ds 

- / F 2 (ph:V h ,T h ,q h ,v h ) + Fi(p h ,v h ,T h , q h ,Th)ds 
Jg 

\ jF> t F 4 (ph, v h , r h , q h , q h ) - 7 \{d t ph)q h j 
Jg 



- F 3 (p hl v h ,T h ,q h ,d t ph)ds 

- / F 2 (ph,v h ,Th,q h ,v h ) + Fi(p h ,v h ,Th,q h ,Th) 
Jg 

+ IPhThVhj ds. 
Thus, an energy conserving scheme has to satisfy 
(3.21) 

0<?J + F a(ph, Vh, Th, q^ dtPh) ds 

,Th) + IPhThVhj ds. 

Note that when (3.211 holds, F4 cannot depend on i>h, t^, q^. Furthermore every 
and F? has to denend on f)tr>i,. As the trace of rttriv. is indenendent 
of the traces of ph, Vh, 'ih, Hh Llie yuaiiuiuies uuiiuaiiimg utp 

and the ones not containing dtPh must cancel each other. 

3.9. Remark. Similarly to the proof of Theorem |3.6| one can show that a scheme 
dissipates energy (even for p = 0), i.e.. 



= / -jD t F 4 (ph,v h ,T h ,q h ,q h ) +jl(d t p h ) 
Jg 

- / F 2 (ph:V h ,T h ,q h ,v h ) + Fi(p h ,v h ,T h ,q h ,' 
Jg 

(3.21| holds, F4 cannot depend on Vh,Th,qh- -furthermore every 

summand in D t F 4 and F3 has to depend on dtPh- As the trace of dtPh is independent 
of the traces of ph,Vh, Th,qh the quantities containing dtPh must cancel each other, 
and the ones not containing dtPh must cancel each other. □ 



(3.22) 



. W ( Ph ) + l - Ph \v h \ 2 + 1 |g,J 2 dx) < 
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provided 



(3.23) 



/ F 1 (p h ,v h ,T h} q h7 T h ) + F 2 {ph,v h ,T h ,q h ,v h ) + lp h T h v h j ds > 
J g 

/ F 3 (p h ,v h ,T h ,q h ,d t ph) - -fB t F i (p h ,v h ,T h ,q h ,q h ) +'jldtp h q h l ds < 0. 
J s 



This opens the way for the construction of schemes which dissipate a small amount 
of energy, which might serve as a stabilising mechanism for the scheme, e.g., in 
case forward time stepping is considered. 

3.10. Example. For a, (3 > the choice of fluxes 

F 1 {p h ,v h ,r h , q h , ¥) = - lp h v h \ {*} +a\r h f I*] 
F 2 (p h ,v h ,r h ,q h ,X) = -^f { p h X} +/? Kl [X] 
F 3 { P h, v h , r h , q h , 3) = -7 {3} 
Fiiph^r^q^Z) = \p h f I Z}, 



(3.24) 



in (3.4 1 yields a scheme which is consistent, conserves mass and dissipates energy. 



Proof of Corollary 3.7 The proof of conservation of mass is exactly the same as in 



the proof of Theorem 3.6 because ( 3.4 1 1 does not depend on p. For the dissipation 



of energy the only difference to the proof of Theorem 3.6 is that, when (3.4)2 is 



tested with X = Vj t an additional summand pB^v^^v^) is created and this term 
is not altered by the subsequent calculations. □ 



We will now show that it is very restrictive for schemes given by (3.4) to be 
consistent, conserve momentum and energy. 

3.11. Proposition. Let be the i-th coordinate vector of \R d . A generic scheme 



d* / PhVh da; = 
In 



of the form (3.4) is momentum conservative, i.e., satisfies 

(3.25) 

if and only if 
(3.26) 

= - / F 1 (p h ,v h ,T h ,q h ,e. l 1 v h ) + F 2 {p h ,v h ,T h ,q h ,e l ) 
J s 

+ F 3 (p h ,v h ,T h ,q h ,d Xi p h ) + F 4 (p h ,v h ,T h ,q h ,Wd Xi p h )ds 

1 



- IPhTheij 



sudn 



W{ Ph ) 

7 IdxiPhQhl + 7 lPh^d XiPh j 

7 



Ph \Vh\ 



- \p h e l 1 (v h <g) Vh) 



1 2 

p h Ap h ei + - \Vph\ e, - d Xi phVph 



ds. 



Proof We start with a rather general calculation of the change of momentum in 



direction ej for i = 1, . . . , d for the generic scheme (3.4). We define 
(3.27) 



Mi{p h ,v h ) := I p h ei J Vhdx 
Jn 



to be the discrete momentum in direction e,- 
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Then the rate of change of the discrete momentum is 



(3.28) 



d t Mi(p h ,v h ) = / d t p h ei J v h + ei J p h dtv h dx. 
Jn 



Now making use of (3.4)i and (3.4| 2 use 4* = e^Vh and X = e», 



(3.29) 

A t Mi(p h ,v h ) = I d\v{p h v h )e i 1 v h - div (p h ej (v h ® v h )) + div (p h v h ) e l 1 v h 
p h ei J WT h + ^p h ei J W (\v h \ 2S ) dx 



/ F 1 (p h ,v h ,T h ,q h ,e i 1 v h ) + F 2 {p hl v h ,T h ,q h ,ei)ds. 
Jg 



Using ( 3.4 13 together with Proposition 3.2 



(3.30) 

d t Mi(p h ,v h ) 



f - div {p h e l 1 (v h <g> v h j) + div (p^e*) r h + ip^eJV (|^| 2 ) da; 

- / Fi(p h ,v h ,T h ,q h ,ei J v h ) + F 2 {ph,v h ,T h ,q h ,ei)ds 
Jg 

- / [Wt^ 



,e, ; 1 ds 



- div {phEi 1 (v h ® «>,)) + div (W(p/,)ej) 

- 7 div (9/,) div (P^e,) + - div (p h \v h | 2 e 4 ) da; 

/ F 1 (p hl v h ,T h ,q h ,e i 1 v h ) + F 2 (p h ,v h .T h ,q h ,ei) 
Jg 

+ F 3 (p h ,v h ,T h ,q h ,div(p fl e i ))ds - / [p^ei] ds. 



Again, in view of Proposition 3.2 we may integrate by parts elementwise and it 
follows 



(3.31) 

d t Mi(p h ,v h ) 



/ IQh^dxtPhdx - / Fi(p h ,v h ,T h ,q h ,e l 1 v h ) + F 2 {p h ,v h ,T h , q h , e t ) 
Jn Jg 



Fz{Ph, Vh, Th, q h , div(p /t e i )) ds 



W^(Ph) + ^Pft l^^l 2 ) e » 



- [p^ej (u/, ® Vh)} - 1 {dxiPhqhl ds - 
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By (3.4)4 and another application of Proposition 3.2 it holds that 



(3.32) 

d t Mi(p h ,v h )= / -^p h Ad Xi p h dx - / F 1 (p h ,v h ,T h ,q h ,e i 1 v h ) + F 2 (p h ,v h ,T h ,q h , 

+ F?,{p h , v h , T h ,q h ,d Xz p h ) + F 4 (p h , v h , T h ,q h ,Wd Xi p h ) ds 

1 



iPhThei 



W( Ph 



:Ph \Vh\ 



ihe-i 1 {vh ® v 



7 IdxiPhQhl + 7 IphVdxtPhl ds. 



We may write the first integral appearing on the right hand side of (3.32) in the 
following way: 

(3.33) 

d t Mi(ph,Vh) = J -7div (^p h Ap h ei + i \V p h \ 2 e 4 - d Xi p h Vp h ^j dx 

- / Fi(Ph,v h ,T h ,q h ,ei J v h ) + F 2 (p h ,v h ,T h ,q h ,ei) 
J s 

+ F 3 (p h , v h ,T h ,q h ,d Xi ph) + F 4 (p h , v h , T h ,q h ,Vd Xi p h ) ds 



is 



{PhTh^il 



sudn 



1 

W(p h ) + -p h \v h \ ) e, 



- IPh&i 1 (v h ® v 



- 7 {dxiPhlJ + 7 {ptNd^phl ds. 
Hence we may apply Proposition |3.2| one more time yielding the desired result. 



□ 



3.12. Remark. Since the right hand side of (3.26) depends nonlinearly on pi ll it 



is not expected that the jump terms involving V ph and q h will necessarily cancel 
with the other terms appearing. Only in the case F4 = will the V ' ph and q h terms 
cancel with each other. In this case (3.9)2 gives us the condition that 



(3-34) F 3 (p h ,v h ,T h ,q h: E) = -7 {q h Sj 

which excludes consistency. In the event that F4 ^ 0, the terms involving Vph in 



(3.26) would be required to cancel independently, yielding a condition on F3 which 



again excludes consistency. 



4. Development of energy consistent numerical methods - the 

temporal discrete case 

For the readers convenience we will present argument for designing the tempo- 
rally discrete scheme in the spatially continuous setting. To obtain a fully dis- 
crete version the spatial and temporal discretisations have to be combined which is 



straightforward as presented in (5.2) and Theorem 5.1 



We subdivide the time interval [0, T] into a partition of N consecutive adjacent 
subintervals whose endpoints are denoted to = < ii < ... < t^ = T. The 
n-th timestep is defined as k n :— t n+ \ — t n . We will consistently use the short- 
hand F n (-) :— F{-,t n ) for a generic time function F. We also denote F n+ z := 
1 (F n + F n+1 ). 
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4.1. Theorem. Given initial conditions p°,v°,T° and q° the temporal semi dis- 
crete scheme is: For n € hi, find p n+1 , v n+1 , t" +1 and q n+1 such that v n+1 = 
and (q n+1 ) J n = on dft and 



(4.1) 



= 



p n+l _ p n 



= p" 



k n 
1 / v n+l 



+ div p n+ 2 v 



div ( p n+ 2 V n+ 2 g t,»+2 ] - div ( p n +2u"+2 ] t,«+2 
2^ 



1 i 1 



- pAv™+2 

1 /! „n+l| 



Kl 2 ) 



= r n + i-^" + ;)-^") +7div ( g n + |)- 
= q«+i _ Vp n+1 . 

This scheme satisfies the following energy dissipation property for all < n < N 



(4.2) 



/ W(p n ) + \p n Kl' + I \q n ? Ax = j n W (P° 



2 

n-l 



1 n I n I 2 7 I n|2 . 

-p |v | + ^ |<Z | dec 



M H fc J / 



Dt; j+ 2 



dec. 



Proof We proceed by multiplying (4.1)! by t 2 and (4.1 ) 2 by v 2, integrate 
over the domain 51 and take the sum. We obtain 



(4.3) 

with 
(4.4) 



0= / f x + f 2 + f + ^4 da; 



fi 



(4.5) 



(4.7) 

f4 



p n+l _ p n _ j^n) ^ / ^ 



p «+l _ 



„n+l 



-p n +2[v n+ 2 



1\T _ w " 



J*2 := div ( p"+2t;"+2 )r" + 2 + p™+2 ( V n ^2 j Vt" + 2 

(4.6) 



1\ T 



f 3 :=(v n+ 2 div p n+ 5»»+5® 1 ,"+ 2 



p™+2 t,"+2 I V 



-div(p n+ 5« n+ t! 
2 



1 

U™+2 



Av n+ 2 dx 
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One may readily check that 



J^ida; 



2 7 



,71+1 I 



(4.8) 



W{p n+1 ) + -p n+1 \v n+L \ + ± \ q 

- JjV{ P n ) + \ P n \v n \ 2 + l\q n \ 2 dx 

- 7 / {p n+1 ~ P n )(q n+1 + q n ) J nds 

Jdn 



dx 



2 7 



W(p n+i ) + -p n+i \v n+L \ + -L \q 



dx 



[ W(p n ) + ~p n \v 



" 9l - n|J r-||q"| 2 dec 



Moreover 
(4.9) 



, 1 , 1 , 1 



^2 dx = div [ p"+2- i ,"+2r" + 2 ) dx= I p n ^2 T n ^2 [ V n ^2 j n ds = 0. 

90 



Furthermore we see that ^ 3 satisfies 



(4.10) 



1 I 1 I 1 \ I 1 



Finally we find for ^4 

?i = P 



(4.11) 



/' 



L 




2 

dx 








L 


T>v n+ l 


2 

dx 









dx-p J (v n+ ^y (r>v n+ ^ nds 



Inserting (|4.8[)-( |4.10| into (|4.3[) yields 
1 
2 

W(p 



= W(p n+L ) + -p 



n+1 |„,n+l| 2 1 T |„i+l| 2 



da; 



(4.12) 
concluding the proof. 



1 



2 p n \v n r + ^\ q >r dx+pk ri 







2 


L 




da;, 









□ 



5. Development of consistent numerical methods - the fully 

discrete case 

In this section we combine our spatial and temporal discretisations to provide 
a fully discrete numerical method for the Euler-Korteweg and Navier-Stokes- 
Korteweg systems. 

Let P v : H 1 ^) -> V,Po : Hj(^) d -> and P~ : H^(^) -> V be the L 2 

o d ri. 

projection operators into V, V and V respectively. We combine the arguments given 
in Sj3] and S|4] to obtain a fully discrete scheme which, given 

(5.1) p° h :=Pyp°, v° h :=P°v°, r°:=P v r and q° h := P~ q° 
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requires us to find a sequence of functions p^ +1 , tJ^ +1 £ V, v'^ 1 g V and <7 ; ™ +1 6 V 
such that 

(5.2) 

= 



+div [p h 2 v h 



^dx 



J $ 



Ph ~ +dlv [Ph v h 



Xdx 



i ,2\T 



Xdx 



( ,. ( n+i n+i\ n+i "+J 1 ri+i 

^-div(p h *e fc 2 )v h -+ Ph -Vr, --^Pft V 

' + * ^^^) +7 div(^)^(K-r + Ki 2 )) s 



pi +i -pi 



da; 



5.1. Theorem. Under the assumptions on the fluxes (3.7), (3.9) given in Theorem 



3.6 the solution of the scheme (5.2) conserves mass, i.e., 



/ pi dx= I pi dx for 0<n<N 



in Jn 
and satisfies the energy dissipation equality 

(5.3) j w(pi +1 ) + \pi +1 K +1 | 2 + 1 1 9 ;; +1 | 2 d* 



l 



Proof The proof is merely combining the results of Theorem 3.6 Corollary 3.7 and 
Theorem ED □ 



6. Numerical experiments 

In this section we conduct a series of numerical experiments aimed at testing 
the robustness of the method. There are four experiments which investigate the 
behaviour of the discrete energy for the Euler-Korteweg (j|6.3[) and the Navier- 



Stokes-Korteweg systems (£6.4), benchmarking the algorithm against a travelling 



wave solution of the Euler-Korteweg system (£6.5), observing that there are no 



parasitic currents in long time simulations. Moreover, we conduct some simulations 



for d = 2 ((6.61 



In each of these experiments we consider the fully discrete scheme (5.2 ) with the 
numerical fluxes given in Corollary |3.8| 
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6.1. Implementation issues. The numerical experiments were conducted using 
the DOLFIN interface for FEniCS [LWlOj . The graphics were generated using 
Gnuplot and ParaView . 

In each of the numerical experiments we fix W to be the following quartic double 
well potential 

(6.1) W(p)= 1 -(p-lf( P -2f 
with minima at p = 1 and p = 2. 

6.2. Remark (the quotient of the double well). In the computational implementa- 



tion we did not use the difference quotient W ^ p n +\_ W n P - appearing in (|4.1|) as it is 



ill-defined for p n+1 = p n and badly conditioned when \p n+1 — p n \ is small. Instead 



we use a sufficiently high order approximation of this term. For (6.1| we use the 
following Taylor expansion representation 

(6.2) W[pn J n 1 + \l^} f>n) = W\p n+ h + ^W"'{p n+ \) - p") 2 

which is exact. We note that when W is not polynomial a sufficiently high order 
truncation of the Taylor expansion can be achieved such that the change in energy is 
of high order with respect to the timestep. This allows the construction of a method 
with arbitrarily small deviations of the energy with respect to the timestep. 

In each of the subsequent numerical experiments we assemble the discrete system 



(5.2) as a nonlinear system of equations. The solution to the nonlinear system was 
approximated by a Newton solver with a tolerance set to 10~ 10 . On each Newton 
step the linear system of equations was approximated using a stabilised conjugate 
gradient solver with an incomplete LU preconditioner also set to a tolerance of 
10 



-10 



6.3. Test 1 conservativity for the Euler Korteweg system. In this case 
we take p = 0. We are then studying the conservativity property of the numerical 
method proposed for the Euler-Korteweg system in Corollary |3.8[ We take f2 = 
[0, 1] and consider an initial condition given by a step function 

(6.3) po{x) = 4 . v = 0. 

I 1.9 otherwise , 

We take 7 = 10 -4 , h — 10~ 4 and k n — k = 10~ 3 for each n. Figure [2] shows the 
energy and mass conservativity of the simulation. 

6.4. Test 2 monotone energy dissipation for the Navier Stokes Korteweg 
system. In this case we take p > and study the dissipation property for the full 



Navier-Stokes-Korte weg system given by (5.2). We take = [0,1] and consider 
the initial conditions (6.3). We fix 7 = 10~ 4 , h — 10~ 4 and k n — k = 10~ 3 . 



We test the effect of the ratio of viscocity to capillarity, i.e., p/"f, on the dynamics 
of the simulation. To that end we run the simulation for p = 10~ 7 (Figure |3|, 
p = IQ- 6 (Figure m and p = 10~ 5 (Figure ^. 
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Figure 2. |6.3| Test 1 - Numerical experiment showing the con- 
servation of mass and energy for the numerical method proposed 



in Corollary 3.8 for the Euler-Korteweg system (i.e., fj, = 0). Due 
to the energy conservativity the Euler-Korteweg simulation will 
never achieve a steady state, the oscillations will continue to pro- 
pogate. 



;i.8 
;i.6 

J. 4 
f 1.2 



1.75 
.1,5 
1.25 



(a) Initial condition, t = 



(b) t = 0.01 




6.5. Test 3 — benchmarking. In this test we look to benchmark the numerical 
algorithm against a steady state solution of the Euler-Korteweg system on the 
domain O = [—1,1]. 
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Figure 3. |6.4| Test 2 - Numerical experiment showing the conser- 
vation of mass and dissipation of energy for the numerical method 
proposed in Corollary |3.8| for the Navier-Stokes-Korteweg system. 
In this test we take /i = 10~ 7 . Notice that the energy dissipation 
allows the Navier-Stokes-Korteweg simulation to achieve a steady 
state. At t — 50 the maximal value of the velocity is of magnitude 
10 -5 . Notice also that /x is chosen sufficiently small such that the 
dynamics are comparible with that of Figure [2] albeit with smeared 
out oscillations. 



;i.e 
;i.6 

J. 4 
f 1.2 



1.75 
1.5 
1.25 



(a) Initial condition, t = 



(b) t = 0.01 




|1.25 
II 

0.996731 




rho 
1.978105 

J .8 

-1.6 

.1.4 

.1.2 

1.021436 



(c) t = 0.05 



(d) t = 0.1 



2.052449 
12 

11.75 



0.1 
0.01 
0.001 
0.0001 ,. 



I otal energy - 
Total mass " 



~K2 CM 0^ 0tB~ 

time 



(e) t = 0.5 



(f) Conservativity plot 



For the double well given by (6.1 ) a steady state solution to the Euler-Korteweg 
system is given by 



(6.4) 
(6.5) 



7 S 3 1 

p(x, t) — tanh , 

FV ' 1 2 2 V2V27 

v(x,t) = Vt 
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Figure 4. |6.4| Test 2 - Numerical experiment showing the effect 
of the ratio of viscocity to capillarity on the dynamics of the simu- 
lation. The simulation is the same as in Figure[3]with the exception 
that /x = 10~ 6 . Notice the oscillations have become smeared out. 
The maximal value of velocity is of mag nitude 1CT 5 at t = 14. 
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(e) t = 0.5 



(f) Conservativity plot 



with appropriate initial data. Note that on the boundary Vp is not zero but 
of negligable value (for small values of 7) . Tables [T]-[3] detail three experiments 
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Figure 5. |6.4| Test 2 - Numerical experiment showing the ef- 
fect of the ratio of viscocity to capillarity on the dynamics of the 
simulation. The simulation is the same as in Figure [3] with the ex- 
ception that /.t = 10~ 5 . Notice the oscillations have become heavily 
reduced in very short time due to the massive dissipation in energy 
initially. The maximal value of velocity is of magnitude 10~ J at 
t = 2. 
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(f) Conservativity plot 



aimed at testing the convergence properties for the scheme for 7 = 10 4 (Tabic [TJ, 
7 = 1(T 5 (Table [2} and 7 = 1(T 6 (Table [3}. 
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Table 1. In this test we benchmark a stationary solution of the 



Euler-Korteweg system using the discretisation (5.2) with piece- 
wise linear elements (p = 1), choosing k — 1/N. This is done by 
formulating ( 5.2 ) as a system of nonlinear equations, the solution to 



this is then approximated by a Newton method with tolerance set 
at 10~ 15 . At each Newton step the solution to the linear system of 
equations is approximated using a stabilised conjugate gradient it- 
erative solver with an successively overrelaxed preconditioner, also 
set at a tolerance of 10~ 15 . We look at the L oo (0, T; L 2 (fi)) er- 
rors of the discrete variables ph and Vh, and use e p := p — ph and 
e v := v — Vh- In this test we choose 7 = 10~ 4 . 



N 




EOC 




EOC 


32 


6.258e-3 


0.000 


6.194e-4 


0.000 


64 


3.028e-4 


4.369 


4.631e-5 


3.742 


128 


4.565e-5 


2.730 


1.105e-5 


2.067 


256 


1.155e-5 


1.983 


3.691e-6 


1.582 


512 


2.945e-6 


1.972 


9.916e-7 


1.896 


1024 


7.368e-7 


2.000 


2.528e-7 


1.972 


2048 


1.842e-7 


2.000 


6.324e-8 


1.999 


4096 


4.605e-8 


2.000 


1.580e-8 


2.009 


2. The test is the same as in Table [T] with the e 


take 7 


= 10~ 5 . 








N 




EOC 


H e "llL 00 (L 2 ) 


EOC 


32 


7.017e-3 


0.000 


1.315e-3 


0.000 


64 


2.469e-3 


1.506 


5.819e-4 


1.176 


128 


4.411e-4 


2.485 


7.672e-5 


2.923 


256 


2.885e-5 


3.935 


5.693e-6 


3.752 


512 


6.5970-6 


2.129 


1.295e-6 


2.136 


1024 


1.668&-6 


1.984 


3.228e-7 


2.004 


2048 


4.1616-7 


2.003 


8.017e-8 


2.010 


4096 


1.0406-7 


2.001 


2.001e-8 


2.003 


3. The test is the same as in Table [T] with the e 


take 7 


= 10~ 6 . 








N 


e P WL2) 


EOC 


e " IIloo(l 2 ) 


EOC 


32 


1.8836-2 


0.000 


1.488e-3 


0.000 


64 


9.071e-3 


1.054 


8.134e-4 


0.871 


128 


3.807e-3 


1.253 


3.820e-4 


1.090 


256 


1.0056-3 


1.922 


9.110e-5 


2.051 


512 


6.486e-5 


3.954 


5.118e-6 


4.171 


1024 


4.907e-6 


3.724 


6.809e-7 


2.910 


2048 


1.016e-6 


2.272 


1.445e-7 


2.236 


4096 


2.439e-7 


2.059 


3.446e-8 


2.068 



26 



JAN GIESSELMANN, CHARALAMBOS MAKRIDAKIS, AND TRISTAN PRYER 



6.6. Test 4 — simulations for d = 2 and parasitic currents. In this test we 
consider the case d — 2. We take S7 = [0, l] 2 and look at the following initial 
condition 



and examine its evolution. 

We expect due to the non-local part of the energy that interfacial layers of size 
~ y/j form, see [Ste88l IORS90 for an energy argument. This process smoothes 
the profile. Moreover, the length of the interface is reduced such that the quadratic 
"droplet" becomes circular. 

We take 7 = p, = 0.0005, h « 0.02 and k n = k = 0.001 for all n. Figure § shows 
the behaviour of the energy and mass of the numerical solution together with the 
solution plot of ph at various times. The solution is overlayed with the velocity Vh 
as a glyph plot. 



(6.6) 
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Figure 6. Test |6.6[ The solution, p h to the Navier-Stokes- 
Korteweg system with initial conditions (6.6 1 at various values of 
t, overlayed with the velocity v^. Notice that the are no parasitic 
currents appearing in the interfacial layer. The velocity tends to 
zero over the entire domain as time increases. The energy-mass 
plot of the simulation is also given. 
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